This example shows how to use R and QGIS from within R to perform a series of common point pattern analysis techniques.
library(mapview)
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(censusxy)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionaddr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/data/master/wic_west_side.csv"))
addr<-addr[c(6, 12:14)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)## street city st zip
## 1 1039 W Hildebrand Ave San Antonio TX 78201
## 2 1102 Delgado St San Antonio TX 78207
## 3 1115 W Martin St San Antonio TX 78207
## 4 1115 W Martin St San Antonio TX 78207
## 5 1115 W Martin St San Antonio TX 78207
## 6 1115 W Martin St San Antonio TX 78207results<-cxy_geocode(addr,
street = "street",
city = "city",
state ="st",
zip = "zip",
class="sf",
output = "simple")## 37 rows removed to create an sf object. These were addresses that the geocoder could not match.results.proj<-st_transform(results,
crs = 2278)OR just use the lat / long information in the data!
addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/data/master/wic_west_side.csv"))
results <- st_as_sf(addr, coords=c("Longitude", "Latitude"), crs=4269,agr="constant")
results.proj<-st_transform(results,
crs = 2278)mv1= mapview(results.proj)
mapshot(mv1, file = paste0(getwd(), "/map1.png"))
library(tmap)
library(tmaptools)
library(OpenStreetMap)
bg<- read_osm(results.proj, ext=1.1)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots()mean feature - average of coordinates
mean_feature<-apply(st_coordinates(results.proj), MARGIN = 2, FUN = mean)
mean_feature<-data.frame(place="meanfeature", x=mean_feature[1], y= mean_feature[2])
mean_feature<-st_as_sf(mean_feature, coords = c("x", "y"), crs= 2278)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green")+
tm_shape(mean_feature)+
tm_dots(col = "red")Central feature - Median of coordinates
median_feature<-apply(st_coordinates(results.proj), MARGIN = 2, FUN = median)
median_feature<-data.frame(place="medianfeature", x=median_feature[1], y= median_feature[2])
median_feature<-st_as_sf(median_feature, coords = c("x", "y"), crs= 2278)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green")+
tm_shape(mean_feature)+
tm_dots(col = "red")+
tm_shape(median_feature)+
tm_dots(col = "blue")Buffer points
wicbuff<- st_buffer(results.proj, dist = 2500)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green")+
tm_shape(wicbuff)+
tm_polygons(col = "red", alpha = .1)Convex hull plot
chull <- st_convex_hull(st_union(results))## although coordinates are longitude/latitude, st_union assumes that they are planartmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green")+
tm_shape(chull)+
tm_polygons(col = "grey", alpha = .5)kernel density - You need projected data for this to work right
R can do kernel density maps, but using simple features it’s kind of complicated. I will use Qgis through R instead using the qgisprocess package
library(qgisprocess)## Using 'qgis_process' at 'C://OSGeo4W64/bin/qgis_process-qgis.bat'.
## Run `qgis_configure()` for details.qgis_configure()## getOption('qgisprocess.path') was not found.## Sys.getenv('R_QGISPROCESS_PATH') was not found.## Trying 'qgis_process' on PATH## Error in rethrow_call(c_processx_exec, command, c(command, args), stdin, : Command 'qgis_process' not found @win/processx.c:994 (processx_exec)## Found 2 QGIS installations containing 'qgis_process':
## C://OSGeo4W64/bin/qgis_process-qgis.bat
## C://OSGeo4W64/bin/qgis_process-qgis-dev.bat## Trying command 'C://OSGeo4W64/bin/qgis_process-qgis.bat'## Success!To use this, we need to find the name of the Qgis algorithm we want. qgis_algorithms() can return all available algorithms, then we can either filter it with View() or use grep to search for one.
algs<-qgis_algorithms()
algs[grepl(pattern = "density", x = algs$algorithm ),]## # A tibble: 8 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:r.li.ed~ r.li.edgedensity r.li.edgedensity
## 2 grass7 GRASS grass7:r.li.ed~ r.li.edgedensit~ r.li.edgedensity.as~
## 3 grass7 GRASS grass7:r.li.pa~ r.li.patchdensi~ r.li.patchdensity
## 4 grass7 GRASS grass7:r.li.pa~ r.li.patchdensi~ r.li.patchdensity.a~
## 5 native QGIS (native c~ native:lineden~ linedensity Line density
## 6 qgis QGIS qgis:heatmapke~ heatmapkernelde~ Heatmap (Kernel Den~
## 7 saga SAGA saga:fragmenta~ fragmentationcl~ Fragmentation class~
## 8 saga SAGA saga:kernelden~ kerneldensityes~ Kernel density esti~qgis_show_help("qgis:heatmapkerneldensityestimation")## C:\Users\ozd504\Documents\GitHub\blog
## Heatmap (Kernel Density Estimation) (qgis:heatmapkerneldensityestimation)
##
## ----------------
## Description
## ----------------
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Point layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## RADIUS: Radius
## Argument type: distance
## Acceptable values:
## - A numeric value
## RADIUS_FIELD: Radius from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## PIXEL_SIZE: Output raster size
## Argument type: number
## Acceptable values:
## - A numeric value
## WEIGHT_FIELD: Weight from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## KERNEL: Kernel shape
## Argument type: enum
## Available values:
## - 0: Quartic
## - 1: Triangular
## - 2: Uniform
## - 3: Triweight
## - 4: Epanechnikov
## Acceptable values:
## - Number of selected option, e.g. '1'
## - Comma separated list of options, e.g. '1,3'
## DECAY: Decay ratio (Triangular kernels only)
## Argument type: number
## Acceptable values:
## - A numeric value
## OUTPUT_VALUE: Output value scaling
## Argument type: enum
## Available values:
## - 0: Raw
## - 1: Scaled
## Acceptable values:
## - Number of selected option, e.g. '1'
## - Comma separated list of options, e.g. '1,3'
## OUTPUT: Heatmap
## Argument type: rasterDestination
## Acceptable values:
## - Path for new raster layer
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT: <outputRaster>
## HeatmapRun the algorithm
wic_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
INPUT=results.proj,
RADIUS = 5000,
PIXEL_SIZE = 100,
KERNEL = 0,
OUTPUT=file.path(getwd(), "wicdenst.TIF"),
load_output = TRUE)## Running "C://OSGeo4W64/bin/qgis_process-qgis.bat" run \
## "qgis:heatmapkerneldensityestimation" \
## "--INPUT=C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN\file215472120b7\file215415475112.gpkg" \
## "--RADIUS=5000" "--PIXEL_SIZE=100" "--KERNEL=0" "--OUTPUT_VALUE=0" \
## "--OUTPUT=C:/Users/ozd504/Documents/GitHub/blog/wicdenst.TIF"
## C:\Users\ozd504\Documents\GitHub\blog
##
## ----------------
## Inputs
## ----------------
##
## INPUT: C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN\file215472120b7\file215415475112.gpkg
## KERNEL: 0
## OUTPUT: C:/Users/ozd504/Documents/GitHub/blog/wicdenst.TIF
## OUTPUT_VALUE: 0
## PIXEL_SIZE: 100
## RADIUS: 5000
##
##
## 0...10...20...30...40...50...60...70...80...90...
## ----------------
## Results
## ----------------
##
## OUTPUT: C:/Users/ozd504/Documents/GitHub/blog/wicdenst.TIFlibrary(raster)
library(RColorBrewer)
result<- qgis_as_raster(wic_dens)
projection(result)<-crs(results.proj)
tmap_mode("view")
tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col="red")+
tm_shape(result)+
tm_raster()Spatial join
A spatial join can combine attributes of one layer with another layer. Here I combine census variables with the WIC clinic points.
library(tidycensus)
library(dplyr)
#load census tract data
sa_acs<-get_acs(geography = "tract",
state="TX",
county = "Bexar",
year = 2019,
variables=c( "DP05_0001E", "DP03_0009P", "DP03_0062E", "DP03_0119PE",
"DP05_0001E","DP02_0009PE","DP02_0008PE","DP02_0040E","DP02_0038E",
"DP02_0066PE","DP02_0067PE","DP02_0080PE","DP02_0092PE",
"DP03_0005PE","DP03_0028PE","DP03_0062E","DP03_0099PE","DP03_0101PE",
"DP03_0119PE","DP04_0046PE","DP05_0072PE","DP05_0073PE",
"DP05_0066PE", "DP05_0072PE", "DP02_0113PE") ,
geometry = T, output = "wide")## Getting data from the 2015-2019 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## Using the ACS Data Profile
## Using the ACS Data Profile#rename variables and filter missing cases
sa_acs2<-sa_acs%>%
mutate(totpop= DP05_0001E, pwhite=DP05_0072PE,
pblack=DP05_0073PE , phisp=DP05_0066PE,
phsormore=DP02_0066PE,punemp=DP03_0009PE, medhhinc=DP03_0062E,
ppov=DP03_0119PE)%>%
dplyr::select(GEOID, totpop, pblack, pwhite, phisp, punemp, medhhinc, ppov)
sa_acs2<-st_transform(sa_acs2, crs = 2278)
sa_trol<-st_cast(sa_acs2, "MULTILINESTRING")spjoin<-st_join(results.proj, sa_acs2)
head(spjoin)## Simple feature collection with 6 features and 229 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 2120097 ymin: 13704380 xmax: 2126362 ymax: 13717570
## projected CRS: NAD83 / Texas South Central (ftUS)
## Source Date Obsolescence.Date
## 1 AtoZDatabases 2/20/2020 8/20/2020
## 2 AtoZDatabases 2/20/2020 8/20/2020
## 3 AtoZDatabases 2/20/2020 8/20/2020
## 4 AtoZDatabases 2/20/2020 8/20/2020
## 5 AtoZDatabases 2/20/2020 8/20/2020
## 6 AtoZDatabases 2/20/2020 8/20/2020
## Business.Name
## 1 Social Workers As Needed
## 2 San Antonio Sign Language Congregation of Jehovah
## 3 Center For Health Care Service
## 4 Linda Hodgson
## 5 Anne Manske
## 6 Adriana T Stanley
## Legal.Name Physical.Address Physical.Address.Number
## 1 1039 W Hildebrand Ave 1039
## 2 1102 Delgado St 1102
## 3 1115 W Martin St 1115
## 4 1115 W Martin St 1115
## 5 1115 W Martin St 1115
## 6 Center For Health Care Service 1115 W Martin St 1115
## Physical.Pre.Direction Physical.Address.Name Physical.Address.Suffix
## 1 W Hildebrand Ave
## 2 Delgado St
## 3 W Martin St
## 4 W Martin St
## 5 W Martin St
## 6 W Martin St
## Physical.Post.Direction Physical.City Physical.State Physical.ZIP
## 1 NA San Antonio TX 78201
## 2 NA San Antonio TX 78207
## 3 NA San Antonio TX 78207
## 4 NA San Antonio TX 78207
## 5 NA San Antonio TX 78207
## 6 NA San Antonio TX 78207
## Physical.ZIP.4 Key.Executive.Name First.Name Middle.Initial Last.Name
## 1 4667
## 2 1506
## 3 3005
## 4 3005 Linda Hodgson Linda Hodgson
## 5 3005 Anne Manske Anne Manske
## 6 3005 Adriana Stanley Adriana Stanley
## Title Gender Location.Employee.Size Corporate.Employee.Size Revenue...Yr
## 1 22 N/A
## 2 17 $1,129,728
## 3 9 $598,091
## 4 Educator Female 1 $66,455
## 5 Educator Female 1 $66,455
## 6 Educator Female NA N/A
## Mailing.Address Mailing.Address.Number Mailing.Pre.Direction
## 1 1039 W Hildebrand Ave 1039 W
## 2 1102 Delgado St 1102
## 3 1115 W Martin St 1115 W
## 4 1115 W Martin St 1115 W
## 5 1115 W Martin St 1115 W
## 6 1115 W Martin St 1115 W
## Mailing.Address.Name Mailing.Address.Suffix Mailing.Post.Direction
## 1 Hildebrand Ave NA
## 2 Delgado St NA
## 3 Martin St NA
## 4 Martin St NA
## 5 Martin St NA
## 6 Martin St NA
## Mailing.City Mailing.State Mailing.ZIP Mailing.ZIP.4 Phone Fax
## 1 San Antonio TX 78201 4667 (210) 521-8100 N/A
## 2 San Antonio TX 78207 1506 (210) 592-6722 N/A
## 3 San Antonio TX 78207 3005 (210) 731-1300 N/A
## 4 San Antonio TX 78207 3005 (210) 358-3731 N/A
## 5 San Antonio TX 78207 3005 (210) 358-3731 N/A
## 6 San Antonio TX 78207 3005 (210) 731-1300 N/A
## Toll.Free County.Name County.Population Metro.Area EIN
## 1 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## 2 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## 3 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## 4 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## 5 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## 6 N/A Bexar 500,000 or more San Antonio-New Braunfels, TX NA
## Main.Line.of.Business Location.Type Importer.or.Exporter
## 1 Health & Allied Services NA
## 2 Deaf Services & Facilities NA
## 3 Counselors & Counseling Services NA
## 4 Counselors & Counseling Services NA
## 5 Counselors & Counseling Services NA
## 6 Counselors & Counseling Services NA
## Manufacturer Primary.SIC Primary.SIC.Description SIC02
## 1 No 8099001 Health & Allied Services 8322011
## 2 No 8322082 Deaf Services & Facilities NA
## 3 No 8322002 Counselors & Counseling Services 8322011
## 4 No 8322002 Counselors & Counseling Services NA
## 5 No 8322002 Counselors & Counseling Services NA
## 6 No 8322002 Counselors & Counseling Services NA
## SIC02.Description SIC03 SIC03.Description SIC04
## 1 Social Workers 9111001 Federal Government Executive Offices 8051005
## 2 NA NA
## 3 Social Workers NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## SIC04.Description SIC05 SIC05.Description SIC06
## 1 Health Care Facilities 8099002 Medical Services Organization 8322087
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## SIC06.Description SIC07 SIC07.Description SIC08 SIC08.Description SIC09
## 1 Clinical Social Workers NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## SIC09.Description SIC10 SIC10.Description NAICS.1
## 1 NA NA NA 621999
## 2 NA NA NA 624190
## 3 NA NA NA 624190
## 4 NA NA NA 624190
## 5 NA NA NA 624190
## 6 NA NA NA 624190
## NAICS.1.Description NAICS.2
## 1 All Other Miscellaneous Ambulatory Health Care Services 624190
## 2 Other Individual and Family Services NA
## 3 Other Individual and Family Services NA
## 4 Other Individual and Family Services NA
## 5 Other Individual and Family Services NA
## 6 Other Individual and Family Services NA
## NAICS.2.Description NAICS.3 NAICS.3.Description NAICS.4
## 1 Other Individual and Family Services 921110 Executive Offices 623110
## 2 NA NA
## 3 NA NA
## 4 NA NA
## 5 NA NA
## 6 NA NA
## NAICS.4.Description NAICS.5
## 1 Nursing Care Facilities (Skilled Nursing Facilities) NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NA
## NAICS.5.Description NAICS.6 NAICS.6.Description NAICS.7 NAICS.7.Description
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## NAICS.8 NAICS.8.Description NAICS.9 NAICS.9.Description NAICS.10
## 1 NA NA NA NA NA
## 2 NA NA NA NA NA
## 3 NA NA NA NA NA
## 4 NA NA NA NA NA
## 5 NA NA NA NA NA
## 6 NA NA NA NA NA
## NAICS.10.Description Non.Profit Number.of.PCs Public...Private Small.Business
## 1 NA No 25 to 99 Private No
## 2 NA No 25 to 99 Private Yes
## 3 NA No 10 to 24 Private Yes
## 4 NA No Up to 4 Private No
## 5 NA No Up to 4 Private No
## 6 NA No Up to 4 Private No
## Square.Footage Website Women.Owned Year.Established
## 1 1,000 to 4,999 www.socialworkersasneeded.com No NA
## 2 1,000 to 4,999 No 2016
## 3 1,000 to 4,999 No NA
## 4 Up to 999 No NA
## 5 Up to 999 No NA
## 6 Up to 999 No NA
## Ticker.Symbol Stock.Exchange Fortune.1000.Ranking Credit.Score
## 1 NA NA NA 0
## 2 NA NA NA 0
## 3 NA NA NA 0
## 4 NA NA NA 0
## 5 NA NA NA 0
## 6 NA NA NA 0
## X2017.Revenue.Yr X2016.Revenue.Yr X2015.Revenue.Yr X2016...Sales.Growth
## 1 0 0 0 0.00
## 2 1129728 1196182 863909 -5.56
## 3 598091 531637 664546 12.50
## 4 0 0 0 0.00
## 5 0 0 0 0.00
## 6 0 0 0 0.00
## X2015...Sales.Growth X2017.Employees X2016.Employees X2015.Employees
## 1 0.00 NA NA NA
## 2 38.46 17 18 13
## 3 -20.00 9 8 10
## 4 0.00 NA NA NA
## 5 0.00 NA NA NA
## 6 0.00 NA NA NA
## X2016...Employee.Growth X2015...Employee.Growth Executive.First.Name.1
## 1 0.00 0.00
## 2 -5.56 38.46
## 3 12.50 -20.00
## 4 0.00 0.00 Linda
## 5 0.00 0.00 Anne
## 6 0.00 0.00 Adriana
## Executive.Last.Name.1 Executive.Title.1 Executive.Gender.1
## 1
## 2
## 3
## 4 Hodgson Educator Female
## 5 Manske Educator Female
## 6 Stanley Educator Female
## Executive.Department.1 Executive.First.Name.2 Executive.Last.Name.2
## 1
## 2
## 3
## 4 Education
## 5 Education
## 6 Education
## Executive.Title.2 Executive.Gender.2 Executive.Department.2
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.First.Name.3 Executive.Last.Name.3 Executive.Title.3
## 1
## 2
## 3
## 4
## 5
## 6
## Executive.Gender.3 Executive.Department.3 Executive.First.Name.4
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.Last.Name.4 Executive.Title.4 Executive.Gender.4
## 1
## 2
## 3
## 4
## 5
## 6
## Executive.Department.4 Executive.First.Name.5 Executive.Last.Name.5
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.Title.5 Executive.Gender.5 Executive.Department.5
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.First.Name.6 Executive.Last.Name.6 Executive.Title.6
## 1
## 2
## 3
## 4
## 5
## 6
## Executive.Gender.6 Executive.Department.6 Executive.First.Name.7
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.Last.Name.7 Executive.Title.7 Executive.Gender.7
## 1
## 2
## 3
## 4
## 5
## 6
## Executive.Department.7 Executive.First.Name.8 Executive.Last.Name.8
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.Title.8 Executive.Gender.8 Executive.Department.8
## 1 <NA>
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## Executive.First.Name.9 Executive.Last.Name.9 Executive.Title.9
## 1
## 2
## 3
## 4
## 5
## 6
## Executive.Gender.9 Executive.Department.9 Executive.First.Name.10
## 1 <NA> NA
## 2 <NA> NA
## 3 <NA> NA
## 4 <NA> NA
## 5 <NA> NA
## 6 <NA> NA
## Executive.Last.Name.10 Executive.Title.10 Executive.Gender.10
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Department.10 Executive.First.Name.11 Executive.Last.Name.11
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Title.11 Executive.Gender.11 Executive.Department.11
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.First.Name.12 Executive.Last.Name.12 Executive.Title.12
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Gender.12 Executive.Department.12 Executive.First.Name.13
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Last.Name.13 Executive.Title.13 Executive.Gender.13
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Department.13 Executive.First.Name.14 Executive.Last.Name.14
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Title.14 Executive.Gender.14 Executive.Department.14
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.First.Name.15 Executive.Last.Name.15 Executive.Title.15
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Gender.15 Executive.Department.15 Executive.First.Name.16
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Last.Name.16 Executive.Title.16 Executive.Gender.16
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Department.16 Executive.First.Name.17 Executive.Last.Name.17
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Title.17 Executive.Gender.17 Executive.Department.17
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.First.Name.18 Executive.Last.Name.18 Executive.Title.18
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Gender.18 Executive.Department.18 Executive.First.Name.19
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Last.Name.19 Executive.Title.19 Executive.Gender.19
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Department.19 Executive.First.Name.20 Executive.Last.Name.20
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Executive.Title.20 Executive.Gender.20 Executive.Department.20
## 1 NA NA NA
## 2 NA NA NA
## 3 NA NA NA
## 4 NA NA NA
## 5 NA NA NA
## 6 NA NA NA
## Est..Accounting.Annual.Expense Est..Advertising.Annual.Expense
## 1 $5,000 to $24,999
## 2 $5,000 to $24,999
## 3 $5,000 to $24,999
## 4 Up to $999
## 5 Up to $999
## 6 Professional Professional
## Est..Business.Insurance.Annual.Expense Est..Legal.Annual.Expense
## 1 $20,000 to $49,999
## 2 $7,500 to $19,999
## 3 $7,500 to $19,999
## 4 Up to $2,499
## 5 Up to $2,499
## 6 Professional Professional
## Est..Office.Equipment.Annual.Expense Est..Rent.Annual.Expense
## 1 $25,000 to $99,999 $25,000 to $99,999
## 2 $10,000 to $24,999 $10,000 to $24,999
## 3 $5,000 to $9,999 $10,000 to $24,999
## 4 Up to $1,999 Up to $2,499
## 5 Up to $1,999 Up to $2,499
## 6 Professional Professional
## Est..Technology.Annual.Expense Est..Telecom.Annual.Expense
## 1 $10,000 to $24,999 $10,000 to $24,999
## 2 $10,000 to $24,999 $5,000 to $9,999
## 3 $5,000 to $9,999 $5,000 to $9,999
## 4 Up to $2,499 Up to $1,999
## 5 Up to $2,499 Up to $1,999
## 6 Professional Professional
## Est..Utilities.Annual.Expense AtoZ.ID Home.Based.Business Franchise.Type
## 1 $20,000 to $49,999 1.11325e+13 No
## 2 $7,500 to $19,999 1.11400e+13 No
## 3 $7,500 to $19,999 1.11325e+13 No
## 4 Up to $2,499 1.11325e+13 No
## 5 Up to $2,499 1.11325e+13 No
## 6 Professional 1.11325e+13 No
## Holding.Parent.Name Source.1 GEOID totpop pblack pwhite phisp
## 1 NA AtoZDatabases 48029190604 4801 1.3 81.8 1.2
## 2 NA AtoZDatabases 48029170101 3881 0.0 87.9 0.7
## 3 NA AtoZDatabases 48029110600 5293 1.1 65.0 6.2
## 4 NA AtoZDatabases 48029110600 5293 1.1 65.0 6.2
## 5 NA AtoZDatabases 48029110600 5293 1.1 65.0 6.2
## 6 NA AtoZDatabases 48029110600 5293 1.1 65.0 6.2
## punemp medhhinc ppov geometry
## 1 1.4 38852 22.9 POINT (2124109 13717566)
## 2 10.4 34375 21.5 POINT (2120097 13707136)
## 3 16.1 15250 37.8 POINT (2126362 13704380)
## 4 16.1 15250 37.8 POINT (2126362 13704380)
## 5 16.1 15250 37.8 POINT (2126362 13704380)
## 6 16.1 15250 37.8 POINT (2126362 13704380)tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(spjoin)+
tm_dots("punemp", size = .1)+
tm_shape(sa_acs2)+
tm_polygons(alpha = .1)Count points in polygons
Point in polygon operations are actually a spatial intersection (more on this next week!) where we see how many points fall within a given polygon.
sa_acs2$nwic<- lengths(st_intersects(sa_acs2, results.proj))
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(sa_acs2)+
tm_polygons("nwic")+
tm_shape(spjoin)+
tm_dots( size = .01)Thiessen/Voronoi Polygons
Thiessen or Voronoi polygons are a process where we can convert points into polygons.
algs[grepl(pattern = "voronoi", x = algs$algorithm ),]## # A tibble: 3 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:v.voronoi v.voronoi v.voronoi
## 2 grass7 GRASS grass7:v.voronoi.sk~ v.voronoi.skele~ v.voronoi.skele~
## 3 qgis QGIS qgis:voronoipolygons voronoipolygons Voronoi polygonsqgis_show_help("qgis:voronoipolygons")## C:\Users\ozd504\Documents\GitHub\blog
## Voronoi polygons (qgis:voronoipolygons)
##
## ----------------
## Description
## ----------------
## This algorithm takes a points layer and generates a polygon layer containing the voronoi polygons corresponding to those input points.
##
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Input layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## BUFFER: Buffer region (% of extent)
## Argument type: number
## Acceptable values:
## - A numeric value
## OUTPUT: Voronoi polygons
## Argument type: sink
## Acceptable values:
## - Path for new vector layer
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT: <outputVector>
## Voronoi polygonswic_von<-qgis_run_algorithm(alg="qgis:voronoipolygons",
INPUT=results.proj,
OUTPUT=file.path(tempdir(), "wicvon.shp"),
load_output = TRUE)
wic_von<-sf::read_sf(qgis_output(wic_von, "OUTPUT"))
tmap_mode("view")
tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(wic_von)+
tm_polygons(col="grey")+
tm_shape(results.proj)+
tm_dots( size = .01)Nearest Neighbor analysis
library(spatstat)
wic.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(wic.pp))
algs[grepl(pattern = "nearest", x = algs$algorithm ),]## # A tibble: 10 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 gdal GDAL gdal:gridinver~ gridinversedist~ Grid (IDW with nea~
## 2 gdal GDAL gdal:gridneare~ gridnearestneig~ Grid (Nearest neig~
## 3 native QGIS (native c~ native:angleto~ angletonearest Align points to fe~
## 4 native QGIS (native c~ native:joinbyn~ joinbynearest Join attributes by~
## 5 native QGIS (native c~ native:nearest~ nearestneighbou~ Nearest neighbour ~
## 6 qgis QGIS qgis:distancet~ distancetoneare~ Distance to neares~
## 7 qgis QGIS qgis:distancet~ distancetoneare~ Distance to neares~
## 8 qgis QGIS qgis:knearestc~ knearestconcave~ Concave hull (k-ne~
## 9 saga SAGA saga:knearestn~ knearestneighbo~ K-nearest neighbou~
## 10 saga SAGA saga:nearestne~ nearestneighbour Nearest neighbourqgis_show_help("native:nearestneighbouranalysis")## C:\Users\ozd504\Documents\GitHub\blog
## Nearest neighbour analysis (native:nearestneighbouranalysis)
##
## ----------------
## Description
## ----------------
## This algorithm performs nearest neighbor analysis for a point layer.
##
## The output describes how the data are distributed (clustered, randomly or distributed).
##
## Output is generated as an HTML file with the computed statistical values.
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Input layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## OUTPUT_HTML_FILE: Nearest neighbour
## Argument type: fileDestination
## Acceptable values:
## - Path for new file
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT_HTML_FILE: <outputHtml>
## Nearest neighbour
## OBSERVED_MD: <outputNumber>
## Observed mean distance
## EXPECTED_MD: <outputNumber>
## Expected mean distance
## NN_INDEX: <outputNumber>
## Nearest neighbour index
## POINT_COUNT: <outputNumber>
## Number of points
## Z_SCORE: <outputNumber>
## Z-scorewic_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
INPUT=results.proj,
OUTPUT_HTML_FILE=file.path(tempdir(), "wicnn.html"),
load_output = TRUE)## Ignoring unknown input 'load_output'## Running "C://OSGeo4W64/bin/qgis_process-qgis.bat" run \
## "native:nearestneighbouranalysis" \
## "--INPUT=C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN\file215472120b7\file21544c7590.gpkg" \
## "--OUTPUT_HTML_FILE=C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN/wicnn.html"
## C:\Users\ozd504\Documents\GitHub\blog## The system cannot find the file C:\rtools40\Version.txt.##
## ----------------
## Inputs
## ----------------
##
## INPUT: C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN\file215472120b7\file21544c7590.gpkg
## OUTPUT_HTML_FILE: C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN/wicnn.html
##
##
## 0...10...20...30...40...50...60...70...80...90...100 - done.
##
## ----------------
## Results
## ----------------
##
## EXPECTED_MD: 1356.5835389412343
## NN_INDEX: 0.4987323168352079
## OBSERVED_MD: 676.5720513566673
## OUTPUT_HTML_FILE: C:\Users\ozd504\AppData\Local\Temp\RtmpawuiPN/wicnn.html
## POINT_COUNT: 100
## Z_SCORE: -9.589602141964955